*** Code for Remotely Incorrect? ****
*** Marginal effects plots

clear all
grstyle init
grstyle set plain


local 	sel_user		2		// 1 - Dan, 2 - Jen
	
	if 		`sel_user' 		== 1 { 
		global 	user		"C:\Users\03638881\Dropbox\Satellite data measurement\"
	}	
	
	if 		`sel_user' 		== 2 {
		global user 		"C:\Users\alixgarj\Dropbox\Satellite data measurement\"
	}
	

	
	global	data			"${user}Data\"
	global 	outputs			"${user}Latex\tables\"
	global 	graphs	 		"${user}Latex\graphs\"
	global 	dofiles			"${user}Do files\"

adopath + "${user}Do Files\"



use "${data}PESi2iPanelData_results_cleanJAG.dta", clear
merge 1:1 yr EjidoCode APEJNCNID using "${data}PESi2iPanelData_results_mc_po_ej.dta"
drop _m
merge 1:1 yr EjidoCode APEJNCNID using "${data}PESi2iPanelData_results_mc_po_poly.dta"
drop _m
save "${data}PESi2iPanelData_results_po_all.dta", replace


***************************
*Plots of marginal effects*
***************************
use "${data}PESi2iPanelData_results_po_all.dta", clear

**Specs 1,2**
preserve
#delimit ;
keep l_1_beneL1 ahl_1_beneL1 s_20_1_beneL1 ahs_20_1_beneL1
l_2_beneL1 ahl_2_beneL1 s_20_2_beneL1 ahs_20_2_beneL1 
mcs_90_11_beneL1 mcl_11_beneL1 mcl_12_beneL1 mcpo_11_beneL1 
mcs_90_12_beneL1 mcpo_12_beneL1 ;
#delimit cr

ren l_1_beneL1 e1
ren ahl_1_beneL1 e2 
ren s_20_1_beneL1 e3 
ren ahs_20_1_beneL1 e4 
ren mcs_90_11_beneL1 e5
ren mcl_11_beneL1 e6
ren mcpo_11_beneL1 e7
ren l_2_beneL1 p1
ren ahl_2_beneL1 p2 
ren s_20_2_beneL1 p3 
ren ahs_20_2_beneL1 p4 
ren mcs_90_12_beneL1 p5
ren mcl_12_beneL1 p6
ren mcpo_12_beneL1 p7

forval i=1/7 {
	qui su e`i', d
	g el`i' = r(p25)
	g eh`i' = r(p75)
	qui su p`i', d
	g pl`i' = r(p25)
	g ph`i' = r(p75)
}
keep if _n==1
keep eh* el* ph* pl* 
g i=1
reshape long eh el ph pl, i(i) j(parm)
replace i=_n
expand 2
bys i: g t=_n
replace eh=. if t==2
replace el=. if t==2
replace ph=. if t==1
replace pl=. if t==1
replace parm=parm -.1 if t==1
replace parm=parm +.1 if t==2
label define parms 1 "CRE Logit" 2 "Ad Hoc CRE Logit" 3 "CRE Scobit" /// 
	4 "Ad Hoc CRE Scobit" ///
	5 "MC-CRE Scobit" 6 "MC-CRE Logit" 7 "PO Probit"
label values parm parms
tw rspike el eh parm, horizontal lc(gs0) lw(vthick) || rspike pl ph parm, horizontal lc(gs8) lw(vthick) xline(0, lc(gs10)) ylab(1(1)7,valuelabel labs(small) angle(0)) ytitle("") xtitle("Marginal Effect") legend(lab(1 "Ejido CREs") lab(2 "Polygon CREs") region(lp(blank))) saving("${graphs}grbox", replace) note("Note: Bars represent the interquartile range of observation-specific marginal effects.", size(vsmall))
graph export "${graphs}grbox_po.png", replace	
restore

loc x1 "everbene elev_met_mean slope_deg_mean dist_any_road_mt dist_5000_city_km TPApAreaHa pctTC2000 per_maj_ind"
glo i=1
foreach x in `x1' {
	preserve
	if ("`x'" == "everbene") loc xx "Ever Beneficiary (0/1)"
	if ("`x'" == "elev_met_mean") loc xx "Average Elevation (meters)"
	if ("`x'" == "slope_deg_mean") loc xx "Average Slope (degree)"
	if ("`x'" == "dist_any_road_mt") loc xx "Distance to any road (meters)"	
	if ("`x'" == "dist_5000_city_km") loc xx "Distance to city with > 5,000 people"	
	if ("`x'" == "TPApAreaHa") loc xx "Area of polygon"
	if ("`x'" == "pctTC2000") loc xx "Percent forested, 2000"
	if ("`x'" == "per_maj_ind") loc xx "Percent of majority indigenous"	
	loc xl = "xline(0, lc(gs10))"
	if ("`x'" == "per_maj_ind" | "`x'" == "pctTC2000") loc xl = " "
	keep l_1_`x' ahl_1_`x' s_20_1_`x' ahs_20_1_`x' mcs_90_11_`x' mcl_11_`x' mcpo_11_`x' 
	ren l_1_`x' p1
	ren ahl_1_`x' p2 
	ren s_20_1_`x' p3 
	ren ahs_20_1_`x' p4 
	ren mcs_90_11_`x' p5
	ren mcl_11_`x' p6
	ren mcpo_11_`x' p7
	forval i=1/7 {
		qui su p`i', d
		g pl`i' = r(p25)
		g ph`i' = r(p75)
	}
	keep if _n==1
	keep ph* pl*
	g i=_n
	reshape long ph pl, i(i) j(parm)
	label define parms 1 "CRE Logit" 2 "Ad Hoc CRE Logit" 3 "CRE Scobit" /// 
		4 "Ad Hoc CRE Scobit" ///
		5 "MC-CRE Scobit" 6 "MC-CRE Logit" 7 "PO Probit"
	label values parm parms
	tw rspike pl ph parm, horizontal lc(gs0) lw(vthick) `xl' ylab(1(1)7,valuelabel labs(small) angle(0)) ytitle("") xtitle("Marginal Effect") title("`xx'", size(*.6)) legend(off) saving("${graphs}grbox${i}", replace) note("Note: Bars represent the interquartile range of observation-specific marginal effects.", size(vsmall))	
	graph export "${graphs}grbox${i}_po.png", replace 
	restore     
	glo i=${i}+1
}


*******************************
*Plots of false negative rates*
*******************************
use "${data}PESi2iPanelData_results_po_all.dta", clear
tw kdensity G0l_11, lc(gs7) lp(solid) || kdensity G1l_11, lc(gs0) lp(solid) || ///
	kdensity G1l_12, lc(maroon) lp(dash) ///
	ytitle("Density", size(*.8)) xtitle("Misclassification Probability", size(*.8)) ylab(,labs(small) angle(0)) ///
	legend(order(- "Ejido FEs:"  1 "  MC-CRE Logit: False Positive" 2 "  MC-CRE Logit: False Negative" ///
	- "Polygon FEs:"  - " " 3 "  MC-CRE Logit: False Negative") col(2) colf size(*.8) ///
	region(lw(none)) rowg(.01) pos(6)) saving("${graphs}gr11", replace) 
graph export "${graphs}mcprates.png", replace	

loc bw=.03
tw kdensity G1s90_11, lc(gs0) lp(longdash_dot) || kdensity G1l_11, lc(maroon) lp(solid) || ///
	kdensity G1po_11, lc(navy) lp(dash) ///
	ytitle("Density") xtitle("Misclassification Probability") ylab(,labs(small) angle(0)) ///
	legend(order(1 "MC-CRE Scobit" 2 "MC-CRE Logit" 3 "PO Probit" ) col(1) colf size(*.8) region(lw(none)) rowg(.01) pos(6)) saving("${graphs}gr12", replace)
graph export "${graphs}mcsrates_po_11.png", replace
loc bw=.03
loc bw1=.08	
tw kdensity G1s90_12, lc(gs0) lp(longdash_dot) || kdensity G1l_12, lc(maroon) lp(solid) || ///
	kdensity G1po_12, lc(navy) lp(dash) ///
	ytitle("Density") xtitle("Misclassification Probability") ylab(,labs(small) angle(0))  ///
	legend(order(1 "MC-CRE Scobit" 2  "MC-CRE Logit" 3 "PO Probit") col(1) colf size(*.8) region(lw(none)) rowg(.01) pos(6)) saving("${graphs}gr13", replace)
graph export "${graphs}mcsrates_po_12.png", replace	

cumul G1l_11, g(cdf_1p)
cumul G1po_11, g(cdf_po)
cumul G1s90_11, g(cdf_90s)
tw connected cdf_90s G1s90_11, sort ms(i) lc(gs0) lp(longdash_dot) || connected cdf_1p G1l_11, sort ms(i) lc(maroon) lp(solid) || ///
	connected cdf_po G1po_11, sort ms(i) lc(navy) lp(dash) ///
	ytitle("Cumulative Probability") xtitle("Misclassification Probability") ylab(,labs(small) angle(0)) yline(1, lc(gs0)) ///
	legend(order(1 "MC-CRE Scobit" 2 "MC-CRE Logit" 3 "PO Probit") colf col(1) size(*.8) region(lw(none)) rowg(.01) pos(6)) saving("${graphs}gr14", replace) 
graph export "${graphs}mcsrates_po_11_cdf.png", replace	
drop cdf_*
cumul G1l_12, g(cdf_1p)
cumul G1po_12, g(cdf_po)
cumul G1s90_12, g(cdf_90s)
tw connected cdf_90s G1s90_12, sort ms(i) lc(gs0) lp(longdash_dot) || connected cdf_1p G1l_12, sort ms(i) lc(maroon) lp(solid) || ///
	connected cdf_po G1po_12, sort ms(i) lc(navy) lp(dash) ///
	ytitle("Cumulative Probability") xtitle("Misclassification Probability") ylab(,labs(small) angle(0)) yline(1, lc(gs0)) ///
	legend(order(1 "MC-CRE Scobit" 2 "MC-CRE Logit" 3 "PO Probit") colf size(*.8) region(lw(none)) rowg(.01) pos(6)) saving("${graphs}gr15", replace) 
graph export "${graphs}mcsrates_po_12_cdf.png", replace	

grc1leg "${graphs}gr12.gph" "${graphs}gr14.gph", row(1) position(6) subtitle("(a) Ejido Fixed Effects")
graph export "${graphs}mcprates_po_11_both.png", replace
grc1leg "${graphs}gr13.gph" "${graphs}gr15.gph", row(1) position(6) subtitle("(b) Polygon Fixed Effects")
graph export "${graphs}mcsrates_po_12_both.png", replace





